#####################################################
# Script written by V. Bartolome on 25 January 2008.#
# Note:  Input file is text file.            Marker #
# names should start with M_ and phenotype names    #
# should start with P_.                             #
#####################################################

inputfile="sample-input-file_osnp20_aus.txt"
outputfile="sample-input-file_osnp20_aus_SMR.txt"
Summary="sample-input-file_osnp20_aus_summary.txt"   # This file can be used to quality check your data.
BarChart="sample-input-file_osnp20_aus.pdf"   # This file shows the distribution of the different marker categories.
alpha=.05  #cutoff for significance, could be tweaked for greater/lesser stringency

data <-read.table(inputfile, sep="\t", header=TRUE)
columnnames=colnames(data)
marker <- grep("M_",columnnames, ignore.case=TRUE, value=FALSE)
pheno <-  grep("P_", columnnames, ignore.case=TRUE, value=FALSE)

mstart=marker[[1]]     
mend=marker[[length(marker)]]
pstart=pheno[[1]]
pend=pheno[[length(pheno)]]

markerdata=data[, mstart:mend]
long=nrow(markerdata)*ncol(markerdata)
markerdata=as.matrix(markerdata)
dim(markerdata)=c(long,1)
markerdata=factor(markerdata)

# set graphics output to pdf
pdf(file = BarChart, paper = "a4"); 
barplot(summary(markerdata), main="Bar Plot for Marker Categories",ylab="Frequency")
dev.off()

for(j in mstart:mend) {
    data[[j]]=factor(data[[j]]) }
capture.output(summary(data), file = Summary)

for(i in pstart:pend) {
   for(j in mstart:mend) {
       level=levels(data[[j]])
       if ((length(level) == 1) && (j==mstart)) mstart=mstart+1
       if (length(level) > 1) {
         model=aov(data[[i]]~data[[j]],data)  ##where ANOVA is done
		 modelTable <- model.tables(model, "means") ## I added, so that the means for each marker level is determined & stored in modelTable object
		 
         phenotype=columnnames[[i]]
         markername=columnnames[[j]]
         F=summary(model)[[1]][["F value"]][[1]]
         P=summary(model)[[1]][["Pr(>F)"]][[1]]
         SSR=summary(model)[[1]][["Sum Sq"]][[1]]
         SSE=summary(model)[[1]][["Sum Sq"]][[2]]
		 ## added by mau - attempt to capture means of M=1, M=3
		 #mark1 <- subset(data[[i]], data[[j]]==1)
		 #mark3 <- subset(data[[i]], data[[j]]==3)
		 m1.mean <-modelTable$tables$'data[[j]]'[["1"]]  ## the mean at marker ==1
		 m3.mean <-modelTable$tables$'data[[j]]'[["3"]] ## mean at marker ==3
		 

         #result=cbind(phenotype,markername,F,P,SSR,SSE)
		 result=cbind(phenotype,markername,F,P,SSR,SSE,m1.mean,m3.mean)  ##includes the means now
         if((i==pstart) && (j==mstart)) overall=result else
         overall=rbind(overall,result)   }
   }
}
overall=data.frame(overall)
overall$phenotype=gsub("P_","",overall$phenotype)
overall$markername=gsub("M_","",overall$markername)
overall$P=substr(overall$P,1,6)
overall$P=as.numeric(overall$P)

overall$F=substr(overall$F,1,6)

overall$SSR=substr(overall$SSR,1,10)
overall$SSR=as.numeric(overall$SSR)

overall$SSE=substr(overall$SSE,1,10)
overall$SSE=as.numeric(overall$SSE)

Rsquare=overall$SSR/(overall$SSE+overall$SSR)
Rsquare=round(Rsquare,6)


overall=cbind(overall,Rsquare)
sig.marker <- overall
sig.marker=subset(overall,P<alpha)
#write.xls(sig.marker, outputfile)
sig.marker$SSR=NULL  ## no need to output
sig.marker$SSE=NULL  ## no need to output
write.table(sig.marker, outputfile, sep="\t", row.names=FALSE, quote=FALSE)

